## This R file Runs the Bayesian Factor Analysis presented in paper
## Code written by Jonathan Slapin
## Final version: October 12, 2011

rm(list=ls())

library(calibrate)
library(foreign)
library(MCMCpack)

# Set working directory
setwd("/Users/jslapin/Dropbox/RTA Exit Options/rio/replication files")
#setwd("/Users/jgray/Documents/My Dropbox/RTA Exit Options/rio/replication files")

# Read in final dataset from Stata
rta<-as.data.frame(read.dta("rta_meanXsal.dta"))

# Subset data to eliminate highly missing RTAs
miss<-ifelse(is.na(rta[,2:ncol(rta)])==T,1,0)
colmiss<-c(0,apply(miss,2,sum))
rta<-rta[,colmiss<20]
rowmiss<-apply(miss,1,sum)
rta<-rta[rowmiss<20,]
rownames(rta)<-rta[,1]
rta<-rta[,-1]

# Bayesian factor analysis
post1dimshort<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"+"),Secretariat=list(1,"+"), IntraPTA=list(1,"+") ), verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=200000,thin=200)


save(post1dimshort,file="post1dim02Aug11.Rdata")

# Summarize Posterior, check for convergence with Geweke diag
sum1dim<-summary(post1dimshort)
gwdiag<-geweke.diag(post1dimshort)
#plot(post1dimshort)


# Calculate percent times RTA phi is greater than EU from posterior sims 
# Keep only phi from posterior

EUphi<-post1dimshort[,40]
postphi1<-post1dimshort[,c(31:39)]
postphi2<-post1dimshort[,c(41:48)]
postphi<-cbind(postphi1,postphi2)

perhighEUbin<-ifelse(EUphi-postphi<0,1,0)
perhighEU<-colSums(perhighEUbin)/1000

# Calculate percent times RTA phi is greater than NAFTA from posterior sims 
# Keep only phi from posterior

NAFTAphi<-post1dimshort[,43]
postphi1<-post1dimshort[,c(31:42)]
postphi2<-post1dimshort[,c(44:48)]
postphi<-cbind(postphi1,postphi2)

perhighNAFTAbin<-ifelse(NAFTAphi-postphi<0,1,0)
perhighNAFTA<-colSums(perhighNAFTAbin)/1000

# Calculate percent times RTA phi is greater than Andean Comm from posterior sims 
# Keep only phi from posterior

Andphi<-post1dimshort[,33]
postphi1<-post1dimshort[,c(31:32)]
postphi2<-post1dimshort[,c(34:48)]
postphi<-cbind(postphi1,postphi2)

perhighAndbin<-ifelse(Andphi-postphi<0,1,0)
perhighAnd<-colSums(perhighAndbin)/1000


##########################

# Extract RTA position scores (phi) from factor analysis
phimeans<-sum1dim$statistics[31:nrow(sum1dim$statistics),1]

# Create Bayesian Credible Intervals
ci<-HPDinterval(post1dimshort,prob=0.95)
philci<-ci[31:nrow(ci),1]
phihci<-ci[31:nrow(ci),2]

# Create dataframe with RTA names and Bayesian Credible Intervals for plotting
rtanames<-cbind("ALADI","APEC","Andean Community","Asean","Bangkok Agreement","CARICOM","CEFTA","ECOWAS","EFTA","EU","GCC","Mercosur","NAFTA","OECS","SACU","SADC","SPARTECA","UEMOA")
bayesciphi<-as.data.frame(cbind(philci,phimeans,phihci))
rownames(bayesciphi)<-rtanames
bayesciphi<-bayesciphi[order(phimeans,decreasing=TRUE),]

# Create plot of 1st dim RTA scores and Bayesian Cred Intervals
pdf(file="firstdim_meansal.pdf")
plot(bayesciphi$phimeans, seq(1:nrow(bayesciphi)),xlab="RTA Positions",ylab="",axes=F,pch=20,xlim=c(-3,2),ylim=rev(range(seq(1:nrow(bayesciphi)))))
axis(1,-1:2)
text(-2.5,seq(1, nrow(bayesciphi)),rownames(bayesciphi),cex=0.8)
segments(bayesciphi$philci, seq(1,nrow(bayesciphi)),bayesciphi$phihci,seq(1,nrow(bayesciphi)))
#abline(v =bayesciphi$phihci[18],lty=2)
dev.off()

# Create matrix of McCall Smith's RTA scores and merge them with our RTA scores
mccsmithrta<-rbind("Andean Community","CARICOM","CEFTA","ECOWAS","EFTA","EU","GCC","Mercosur","NAFTA","OECS","SACU")
mccsmithscore<-t(as.matrix(cbind(4,1,0,3,0,4,1,2,2,2,0)))
rownames(mccsmithscore)<-mccsmithrta
ndf<-merge(mccsmithscore,bayesciphi, by.x=0,by.y=0)
ndf$Row.names[1]<-"Andean Com"


# Plot McCall Smith against our dimension
pdf(file="plotMcCall.pdf")
fit<-lm(ndf$phimeans~ndf$V1)
cor(ndf$V1,ndf$phimeans)
plot(ndf$V1,ndf$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1),xlab="McCall Smith Legalization Score",ylab="Expert Survey Principal Dimension")
textxy(ndf$V1,ndf$phimeans,ndf$Row.names,cx=.75)
abline(fit)
dev.off()

# Plot BearceOmori and Grigorescu dimensions against ours.
ndf2<-read.dta("comparisons.dta")
ndf2$phimeans<-bayesciphi$phimeans
fit<-lm(phimeans~integration,data=ndf2)
fit2<-lm(phimeans~integration,data=ndf2[-10,])
cor(ndf2$integration,ndf2$phimeans)

pdf("BearceOmori.pdf")
plot(ndf2$integration,ndf2$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1),xlab="Bearce-Omori Integration Score",ylab="Expert Survey Principal Dimension")
textxy(ndf2$integration,ndf2$phimeans,ndf2$rtaname,cx=.75)
abline(fit)
abline(fit2,lty=2)
dev.off()

fit<-lm(phimeans~oversightfunctions,data=ndf2)
fit2<-lm(phimeans~ oversightfunctions,data=ndf2[-10,])
cor(ndf2$oversightfunctions,ndf2$phimeans)
ndf3<-ndf2[is.na(ndf2$oversightfunctions)!=T,]

pdf("Grigorescu.pdf")
plot(ndf3$oversightfunctions,ndf3$phimeans,pch=20,xlim=c(0,4.8),ylim=c(-.55,1.1), xlab="Grigorescu Transparency Score",ylab="Expert Survey Principal Dimension")
textxy(ndf3$oversightfunctions,ndf3$phimeans,ndf3$rtaname,cx=0.75)
abline(fit)
dev.off()

# Run different chains for Bayesian Factor analysis with different starting values to check for convergence

post1dimshortc2<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"+")), lambda.start=.3,psi.start=-.5,verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=150000,thin=150)

post1dimshortc3<- MCMCfactanal(~Ability + Ambition + Delegation + DSM + Enforcement + Goods + IntraPTA + Legalization + Monitoring +  NTBs  + Scope + Secretariat + Services +  Tariffinternal + Tradediversion,factors=1,lambda.constraints=list(Ability=list(1,"-")), lambda.start=runif(1,-.25,.25), psi.start=runif(1,-.25,.25),verbose=0, store.scores=TRUE, a0=1, b0=0.15, data=rta, burnin=50000,mcmc=150000,thin=150)

